Hierarchy of Fast Motions in Protein Dynamics 
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For many biological applications of molecular dynamics (MD) the importance of good sampling in 
conformational space makes it necessary to eliminate the fastest motions from the system in order to 
increase the time step. An accurate knowledge of these motions is a necessary prerequisite for such 
efforts. It is known that harmonic vibrations of bond lengths and bond angles produce the highest 
frequencies in proteins. There are also fast anharmonic motions, such as inter-atomic collisions, 
which are probably most important when bond lengths and bond angles are fixed. However, the 
specific time scales corresponding to all these limitations are not known precisely. In order to clarify 
the above issue this paper analyses time step limiting factors in a series of numerical tests by using 
an internal coordinate molecular dynamics approach, which allows chosen internal coordinates to 
be frozen. It is found that, in proteins, there is a rather complicated hierarchy of fast motions, 
with both harmonic and anharmonic effects mixed together at several close time scales. Non- 
bonded interactions, notably strong hydrogen bonding, create locally distributed normal modes 
with frequencies similar to those of bond stretching between non-hydrogen atoms. They also impose 
ubiquitous anharmonic limitations starting from rather small step sizes. With fixed standard amino 
acid geometry, rotations of hydrogen bonded hydroxyl groups limit time steps at the 5 fsec level. The 
next important limitation occurs around 10 fsec and is created by collisions between non-hydrogen 
atoms. 



I. INTRODUCTION 

In the literature devoted to protein dynamics two 
classes of motions are most frequently discussed. The 
first comprises the slowest motions involved in the fold- 
ing and various biological functions of proteins. The 
second class is at the opposite border of the spectrum 
and includes the fastest motions which limit time steps 
in molecular dynamics (MD) calculations and, conse- 
quently, the time scales accessible for simulation. Be- 
cause of the high computer cost of long MD trajectories 
for biological systems, considerable efforts have always 
been applied to overcome these linutatio-ns.-.The first and 
most popular method is constraint MDErati. Among the 
techniques developed in recent years-one can mention, 
for example, multiple time scale MDoEl, weig 
MDEH3 and internal coordinate MD (ICMD) 
Within the framework of this research the physical na- 
ture and the specific time scales aLthe fast motions in 
proteins have often been discussecfj. There are several 
reasons, however, why this issue required the more com- 
plete analysis presented here. 

One of the motivations came from recent progress in 
ICMD, that is dynamics simulations in thcpspacc of in- 
ternal rather than Cartesian coordinatesEJllj. By this 
method, proteins can be modeled as chains of rigid bod- 
ies with only torsional degrees of freedom. It is not ob- 
vious a priori what are the time step limiting motions 
in such models, and, in the literature, there are rather 
contradictory estimates of the-mffirall prospects of ICMD 
compared to constraint MDEflO. Another reason why 
this question attracts attention came from the recent 



results concerning the properties of the numerical algo- 
rithms commonly employed in MD-jaamcly integrators of 
the Stormer-Verlet-leapfrog groupli3. It has been found 
that the fluctuation of the computed instantaneous to- 
tal energy of a microcanonical ensemble, which is gener- 
ally used as the most important test for quality of MD 
trajectories, is dominated by simple interpolation errors 
and cannot serve as a measure of accuracy. It is possible, 
therefore, that previous estimates of limiting step sizes 
were significantly biased. 

Proteins have hierarchical spatial organization and this 
structural hierarchy is naturally mapped onto the spec- 
trum of their motions. Namely, fast motions involve in- 
dividual atoms and chemical groups while the slow ones 
correspond to displacements of secondary structures, do- 
mains etc. Every such movement considered separately 
can be characterized by a certain maximum time step, 
and in this sense one can say that there exists a hierar- 
chy of fast motions and, accordingly, of step size limits. 
The lowest such limit is determined by bond stretching 
vibrations of hydrogens, but for our purposes the fol- 
lowing few levels of this hierarchy are most interesting. 
Normally it is assumed that stretching of hands between 
non-hydrogen atoms forms the next leveltj. However, 
this intuitive suggestion is not easy to verify because pro- 
tein normal modes in this frequency range are not always 
local, and their high frequencies cannot be readily at- 
tributed to specific harmonic terms in the force field. On 
the other hand, normal modes with bond lengths con- 
strained have never been studied. It should also be ex- 
pected that, because of the limitations imposed by fast 
anharmonic motions, the hierarchy does not exactly cor- 
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respond to the spectra of normal modes. All these issues 
are considered in detail below. 

The paper is organized as a simple sequence of numer- 
ical tests on a model protein with the fastest motions 
suppressed one after another. Step size limits are de- 
termined by a-unique test proposed and analyzed in de- 
tail elsewhereta. It is found that the hierarchy of step 
size limits is rather complicated and does not always fol- 
low common intuitive assumptions. For instance, bond 
stretching between non-hydrogen atoms in fact overlaps 
the frequency range of certain collective vibrations and, 
therefore, does not create a separate limitation on time 
steps. On the other hand, with step sizes beyond 10 fsec, 
anharmonic effects become dominant. In agreement with 
recent studies£3, leapfrog trajectories appear to hold to 
correct constant-energy hypersurfaces with considerably 
larger time steps than commonly recommended. 

II. RESULTS AND DISCUSSION 
Test Systems and Simulation Protocols 

The main model system consists of an immunoglobulin 
binding domain of streptococcal protein Go which is a 
56 residup-p//3 protein subunit (file lpgb in the Protein 
DatabaseEl) with 24 bound water molecules available in 
the crystal structure. All hydrogens are considered ex- 
plicitly giving a total of 927 atoms. In order to separate 
effects produced by water from those due to the protein 
itself similar tests were also performed for a droplet of 
TIP3P water molecules. The initial configuration of the 
droplet was obtained by taking the coordinates of the 
first 100 molecules closest to the center of the,|Standard 
water box supplied with the AMBER packagec3E3. Equi- 
libration for low temperature tests was performed start- 
ing from the local energy minimum corresponding to this 
configuration rather than from an ice crystal structure. 

In all ^ealculations the AMBER94 force field was 
employedc-3 without truncation of non-bonded interac- 
tions. Molecular motions were generally simulated with 
internal coordinates as independent variables by using 
Hamiltonian equations of motios. and an implicit leapfrog 
integrator described elsewherdla In one case, however, 
namely for normal temperature calculations with fixed 
bond lengths to hydrogen atoms, the standard AMBER 
package was employed, because the ICMD algorithm 
needs too_ many iterations for convergence with large 
step sizesEj. Comparisons between ICMD trajectories 
and usual Cartesian MD show no essential differencesllj, 
therefore, for consistency and convenience, internal coor- 
dinates have been preferred wherever possible. 

Initial data for all numerical tests wece-prepared with 
the standard protocol described earlieni§ta which makes 
possible smooth initialization of leapfrog trajectories al- 
ways from a single constant energy hypersurface. In or- 
der to simulate a virtually harmonic behavior parallel 



tests were performed at the very low temperature of IK 
for which the equilibration protocol was modified as fol- 
lows. During the first 5 psec all velocities were reas- 
signed several times by sampling from a Maxwell distri- 
bution with T=l K. During the following 7.5 psec ve- 
locities were rescaled periodically if the average temper- 
ature went above 2 K. This modification was necessary 
since in the virtually harmonic low temperature condi- 
tions energy equipartition is reached very slowly and the 
initial distribution over normal modes can persist for a 
long time. The necessary harmonic frequencies were ob- 
tained from calculated spectral densities of autocorrela- 
tion functions of appropriate generalized velocities. In all 
production runs the duration of the test trajectory was 
10 psec. 

The estimates of maximal time steps are normally 
made by checking consecration of the total energy of a mi- 
crocanonical ensembleEjcJ. Our approach is similar, but 
it accurately takes into account certain non-trivial prop- 
erties of the leapfrog discretization of MD trajectorieala. 
The test trajectory is repeatedly calculated always start- 
ing from the same constant-energy hypersurface. In each 
run certain system averages are evaluated and compared 
with "ideal" values, i.e. the same parameters obtained 
with a very small time step. The choice of such pa- 
rameters must correspond to the leapfrog discretization, 
which means that they should be computed from on-step 
coordinates and half-step velocities without additional 
interpolations. The latter condition, together with the 
"smooth start" , distinguishes this approach from ear- 
lier testing strategies. These modifications are essential 
because they remove a significant and systematic bias 
present in the traditional approach which employs the 
time fluctuation of thcj instantaneous total energy as an 
indicator of accuracytj. 

The parameters we use are as follows: the average po- 
tential energy, U, and its time variance, D[U]; the aver- 
age kinetic energy, K, computed for half-steps; the total 
energy, E = U + K , and its drift computed in the same 
way, referred to below as _E-drift. For a sufficiently long 
trajectory of a Hamiltonian system U, E and D[U] char- 
acterize the sampled hypersurface in phase space. Their 
deviations from the corresponding virtually ideal values 
characterize the bias of the sampling obtained and, there- 
fore, can be used for accessing step size limits. The E- 
drift computed-4n this way is exactly zero for ideal har- 
monic systems^, and is thus a good indicator of anhar- 
monic effects. 

As an example let us consider results of such testing 
for a completely free protein. It is seen in Fig. [l] that, 
below a step size of about 1.7 fsec, all the measured pa- 
rameters remain approximately constant and close to the 
accurate values. Above this level some deviations grow 
rapidly. This characteristic behavior is s im ilar to that of 
a simple leapfrog harmonic oscillatorli§E3. All its prop- 
erties depend upon the reduced step size r = toh, where 
u> is the frequency, and it can be shown analytically that, 
with small r, power series expansions for deviations of E 
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and Uj&re dominated by the terms of the fourth and sixth 
orderscS. This explains why the deviations grow rapidly 
beyond a certain threshold which depends mainly upon 
the highest frequencies in the system. 
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FIG. 1. Characteristic time step dependencies obtained for 
a completely free protein. The corresponding low temperature 
plots are shown by thinner lines in (a-c). For comparison, in 
(a) and (c) the low temperature deviations for the total and 
potential energies had been scaled and shifted to fit the range 
of deviations observed with normal temperature. The dashed 
horizontal lines in (a) and (c) show the bands of acceptable 
deviation corresponding to variances D[U] indicated by the 
dashed line in (d). The definition of these bands is given in 
the text. 

In order to access the accuracy quantitatively we need 
to compare the deviations of average energies with some 
scale. An appropriate natural scale is given by the time 
variations of the potential energy characterized by the 
value of D[U] shown in Fig. [I] (d). We will take a devi- 
ation of 0.1D[i7] as the upper acceptable level for U and 
a two times larger value for E. In a harmonic system, U 
and K computed as described are equally ; so the devia- 
tion of E is exactly two times that of U and it reaches 
its upper level simultaneously with U , i.e. with the same 
characteristic step size denoted as h c . The threshold lev- 
els chosen are certainly arbitrary, but they are reasonable 
and normally h c values appear similar to maximal step 
sizes reported in the literature. We note, however, that 
we will be mainly interested here in relative rather than 
in absolute h c values for different models. 

Thus, the dotted line in Fig. [j] (d) marks the best 
estimate of the variance D[U] and similar lines in Figs. 
|l| (a,c) show the corresponding acceptance intervals for 
energies. It can be seen in Figs. [I] (a) and (c) that, at 
normal temperature, the deviations of the total and po- 
tential energies look qualitatively different because U is 
affected by occasional transitions between local minima 



which are stochastic and cannot be properly averaged 
during the relatively short test trajectory. Nevertheless, 
both U and E go beyond their acceptance intervals with 
h c rs 1.7 fsec. At low temperature, both deviations are 
regular, and they yield almost exactly the same h c value. 
The highest frequencies in this system are those of the 
bond stretching modes of hydrogens and they range from 
approximately 3000 cm' 1 for aliphatic-scoups to 3800 
cm -1 for the fastest hydroxyl groups£irc3c3. It appears, 
therefore, that our h c value corresponds to t « 1.1 pLn 
reasonable agreement with an ideal harmonic modclEa. 
For a frequency of 3800 cm -1 the stability limit of the 
leapfrog scheme is h = 2.78 fsectll, and in low tempera- 
ture tests the numerical trajectory remains stable right 
up to this value. At normal temperature, however, a sig- 
nificant E-drift appears with h > 2 fsec, which indicates 
that some fast anharmonic motions occur in the system, 
and the test trajectory could not be completed because 
of an explosive growth in temperature. The high E-drift 
can be caused, for instance, by collisions between hydro- 
gens in non-polar contacts. 

Since the data presented in Fig. [l] are redundant each 
case below is characterized by plots (a) and (b) alone. 
The deviation of the total energy is sufficient to evaluate 
h c , while U always behaves similarly to plot (c). As for 
the value of D[U] shown in Fig. [j] (d), its large and 
systematic deviation always results from the drift of the 
potential energy which is implicitly included in plot (b). 
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FIG. 2. Protein models with bond length constraints: (a,b) 
- only bonds to hydrogen atoms fixed; (c,d) - all bond lengths 
fixed. The notation is similar to Figs, hi (a) and (b). 



Proteins with Constrained Bond Lengths 

The two standard modes of constraining bond lengths, 
that is, constraining only bonds to hydrogens and con- 
straining all bonds, yield the results shown in Fig. 
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We see that they exhibit rather similar behavior with h c 
values around 3 fsec for both low and normal tempera- 
ture simulations. In proteins, the fastest bond stretching 
modes between non- hydrogen atoms normally occur jn 
carboxyl groups, with frequencies around 1720 cm -1 E3. 
In our calculations, however, the fastest mode was found 
at 1850 cm -1 a tryptophan side chain, which, in a har- 
monic case, would give a maximum step size of 5.7 fsec 
and h c w 3. We note, therefore, that the plots in Figs. |^ 
(a,b) agree well with a harmonic approximation and give 
the expected h c value with only bonds to hydrogens con- 
strained. With all bonds constrained, however, the ori- 
gin of the limitation is not clear. The fastest bond angle 
bending modes are around 1600 cm^ 1 cB, i.e. not so far, 
but still they should not impose limitations in this range 
of time steps. These limitations are evidently harmonic, 
however, because in low temperature tests (see Figs, g 
(a) and (c)) the two models behave almost identically. 
Note that they both almost reach the harmonic step size 
limit, with E-drift close to zero, while the h c values are 
roughly the same as those at normal temperature. 
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FIG. 3. Spectral densities of velocity autocorrelation func- 
tions for H-N e -H bond angle of Lysl3 for low tempera- 
ture MD trajectories calculated with only bonds to hydrogen 
atoms fixed (a), and with all bond lengths fixed (b). The 
spectra had been computed from separate 40 psec runs with 
1 fsec time step and all points stored. The densities are given 
in arbitrary scales. 

An explanation to this apparently counter-intuitive be- 
havior follows from the example shown in Fig. ||. This 
figure presents spectral densities of bond bending oscil- 
lations of one of the H-N e -H angles of Lysl3 obtained 
at low temperature with the two different modes of ap- 
plying bond length constraints. We see that, with fixed 
bonds to hydrogen atoms only, the maximum occurs at 
1780 cm -1 , not very far from the fastest bond stretch- 



ing. These two hydrogens are involved in hydrogen bonds 
with the peptide oxygen of Gly9 and a water molecule, 
respectively. Neither this angle, nor any of the adjacent 
bonds or bond angles have independent frequencies above 
1700 cm -1 , nevertheless, the signal at 1780 cm -1 is ob- 
served in the spectral densities of many valence and dihe- 
dral angles between neighboring atoms, as well as within 
the hydrogen bonded Gly9, indicating that this frequency 
corresponds to a locally distributed normal mode. 

With all bond lengths fixed, this normal mode can- 
not remain intact, but it does not disappear, which is 
clear from Fig. || (b). Compared with Fig. || (a) there 
are many fewer signals above 1500 cm -1 , but one at 
1690 cm^ 1 has appeared. This peak apparently corre- 
sponds to the same fast collective mode as in Fig. ^ (a), 
with bond stretching eliminated, and, as we see, its fre- 
quency is just slightly reduced. Such behavior is rather 
characteristic of bond angle vibrations of hydrogens in- 
volved in hydrogen bonding, which explains similar step 
size limitations for the two modes of bond length con- 
straint. A high frequency in this case results from a 
combination of several terms in the force field, rather 
from a single specific one, with hydrogen bonding as one 
of the major components. Bond length constraints affect 
such modes indirectly, mainly due to redistribution the 
system inertia. We note finally that this particular exam- 
ple had been selected because in both constraint modes 
there is only one peak above 1650 cm -1 , therefore, in the 
spectrum shown in Fig. ^ (a), it cannot be attributed to 
bond stretching. The largest independent shifts due to 
hydrogen bonding are observed in valence angle vibra- 
tions of hydroxyl groups and their frequencies also reach 
the level of 1700 cm -1 . 



Water Droplet 

The results of tests with a droplet of TIP3P water 
molecules are shown in Figs. |] (a,b). They exhibit an 
h c value of 5 fsec, although one should note that beyond 
6 fsec, with normal temperature, the E-drift grows very 
rapidly. It might be expected that, since in this case no 
purely harmonic terms are present in the force filed, the 
limiting motion should be anharmonic. The low temper- 
ature tests, however, yield exactly the same h c value and, 
thus, it appears that the system behaves rather similarly 
to the previous ones with a time step limiting frequency 
of approximately 1100 cm -1 . This value is close to the 
experimental upper boundary of the band .attributed to 
the rotations of individual water moleculesEj. 

Rotation, or rather libration, of a single water molecule 
in a net of hydrogen bonds is certainly the fastest motion 
here and it can be specifically slowed down by artificially 
increasing the moments of inertia of the water molecule. 
Such water models present considerable interest for sim- 
ulations where structural and thermodynamic properties 
are targeted, rather than kinetic onesa. Figures El (c,d) 
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demonstrate results of such tests with an inertia Iy = /i(5ij 
added to oxygen atoms, where Sij is Kronecker delta and 
fx = 4 (atom mass units) -A 2 . This means that oxygens 
are no longer considered as point masses but as spherical 
rigid bodies of the same mass. With // = 4 the highest 
frequency is expected to be reduced approximately by a 
factor of two, and it is seen that plots in Figs. (c,d) 
look like those in Figs. || (a,b) scaled along the hori- 
zontal axis, leading to a twofold increase in h c . Further 
increase of inertia gives the effect shown in Fig. |J (e,f) 
where [i equals 15. Instead of an increase proportional to 
the square root of the added inertia we obtain h c around 
14 fsec and 10 fsec for low and normal temperatures, 
respectively. In this case, therefore, we encounter a qual- 
itatively different situation with significantly anharmonic 
limitations probably imposed by the translational motion 
of water molecules and collisions between them. 
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FIG. 4. Water droplet tests. Results are shown for three 
different water models: (a,b) - standard TIP3P water model; 
(c,d) - added inertia equaled 4; (e,f) added inertia equaled 15. 
The notation is as in Figs. ^ (a,b). 

By considering the results shown in Fig. |J one con- 
cludes that the water model with fj, = 4 is "well bal- 
anced" for normal temperature simulations in a sense 
that both its translational and rotational movements oc- 
cur in the same time scale and require similar time steps. 
It is worth noting that the step size limits obtained here 
with the new test agree well with earlier results. It is 
known, for instance, that, with a rigid water model, all 
liquid structural properties are accurately reproduced up 
to a step size of 6 fsectiil. A similar assertion holds for 
the "weighted mass" water model up to a step size of 
10 fsecfl. 



Proteins with Fixed Standard Amino Acid 
Geometry 

In these calculations bond lengths and bond angles 
in the protein were fixed according to a standard ge- 
ometry approximations. The results shown in Figs. 
^ (a,b) were obtained with no modifications of inertia 
tensors. It is seen that they strongly resemble the water 
droplet plots in Figs. [| (a,b). This might mean that the 
step size limitations are imposed by a few water clusters 
around the protein, but calculations with increased wa- 
ter inertia yield virtually identical results. On the other 
hand, simultaneous weighting of water and protein hy- 
droxyl groups yields a considerable effect as shown in 
Figs. [| (c,d). The additional inertia tensor was same as 
in the water droplet tests above. Thus, it is evident that, 
in the standard geometry representation, libration of hy- 
droxyl groups is the fastest motion, certainly due to their 
small inertia. The high frequency of these librations is 
due to hydrogen bonding rather than to the correspond- 
ing torsional potential which produces oscillations with 
at least three times lower frequencies. 
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FIG. 5. Calculations with constraints corresponding to the 
standard amino acid geometry. (a,b) - with no modifications 
of inertia tensors; (c,d) - with inertia of water molecules and 
protein hydroxyl groups increased by 4; (e,f) - with inertia 
tensors of water molecules and protein hydroxyl and amide 
groups increased by 36, and masses of water oxygens equal to 
72. The notation is as in Figs, [i] (a,b). 

We saw above that a weighted water model with fi = 4 
gives h c « 10 fsec, and that this value already does not 
depend upon the dynamics of hydrogen motion. It seems 
that, in general, the 10 fsec level is characteristic of sys- 
tems in which the time scale difference between the move- 
ments of hydrogen and non-hydrogen atoms is somehow 
smoothed, which is illustrated by Figs. || (e,f). In these 
calculations, in order to remove all the limitations im- 
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posed by water, the masses and inertia tensors of wa- 
ter oxygens were increased by 56 and 36, respectively, 
which gives a fourfold increase of both masses and iner- 
tia tensors compared with Figs. [| (c,d) and, accordingly, 
a twofold increase in step sizes. The same inertia was 
added to protein hydroxyl and amide groups. In the low 
temperature plots in Figs. || (e,f), it is seen that the har- 
monic h c is shifted compared to Fig. |^ (c), which means 
that, in the previous model, librations of hydroxyl and 
amide groups produce the highest frequencies. The effect, 
however, takes place at low temperature only, and the 
characteristic time step in Fig. || (e) is approximately the 
same as in Fig. ^| (c). A similar behavior is observed in 
calculations with weighted inertia of other hydrogen-only 
rigid bodies (results not shown). The situation, there- 
fore, appears very similar to that in the water droplet 
tests in Fig. |J, but in the present case fast collisions be- 
tween non-hydrogen protein atoms present the limiting 
factor. A possiblq-jway to overcome this limitation is the 
RESPA approachO which until now has been used only 
within the context of Cartesian MD, but cajj-.be adapted 
for ICMD with the Hamiltonian equations^] since they 
make possible symplectic numerical integration. 

III. CONCLUDING REMARKS 

Although the problem of time step limitations in MD 
simulations of biopolymcrs is long standing, and although 
the search for effective remedies has been the subject of 
many studies, a precise account of the effects that can 
create such limitations was missing in the literature. It 
is necessary to fill this gap in order that efforts put into 
the development of new methods can be effective. 

The conventional view of this problem is that the step 
size is limited by fast harmonic motions and that these 
motions are produced by the stiff components of empiri- 
cal potentials connected with the deformations of bonds, 
bond angles, planar groups etc. It has been shown here 
that this description takes into account only a part of step 
limiting factors. First, there may be many fast harmonic 
vibrations in which hydrogen bonding plays a major role. 
Second, non-bonded inter-atom interactions impose ubiq- 
uitous anharmonic limitations starting from rather small 
step sizes. All such factors form a complicated hierarchy 
which slightly differs between different empirical poten- 
tials and which is modified when constraints are imposed. 

Our results suggest that the time steps currently used 
in MD simulations can be increased considerably with no 
serious loss in the accuracy of thermodynamic averages. 
In unconstrained calculations with AMBER potentials 
the estimated limit is about 1.7 fsec. Note, however, 
that this value is mainly connected with the frequency 
of OH-stretching which is high in AMBER because it 
close to the upper limit of experimental frequencies of 
free hydroxyl groupsEa. In the ENCAD potentials^, for 
instance, the corresponding force constant is lowered by 



a factor of 2, which is why a larger step size of 2 fsec 
recommended by the authors is safe and could possibly 
even be increased further. 

Constraining bonds to hydrogen atoms removes the 
lowest level in the hierarchy of fast motions and makes 
possible time steps up to 3 fsec. Constraints on other 
bond lengths, however, are effective only together with 
constraining bond angles. With fixed standard amino 
acid geometry, rotations of hydrogen bonded hydroxyl 
groups limit time steps at a 5 fsec level. The next im- 
portant limitation occurs around 10 fsec and it is due to 
collisions between non-hydrogen protein atoms. 
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